{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "aIxf_rKkFDsg"
   },
   "source": [
    "<center>\n",
    "<img src=\"https://habrastorage.org/files/fd4/502/43d/fd450243dd604b81b9713213a247aa20.jpg\" />\n",
    "    \n",
    "## [mlcourse.ai](mlcourse.ai) – Open Machine Learning Course \n",
    "### <center> Author: Syrovatskiy Ilya Igorevich, ODS Slack nickname : bokomaru\n",
    "    \n",
    "## <center> Individual data analysis project \"Rossmann store sales\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "N283Sr6VFDsh"
   },
   "source": [
    "**Research plan**\n",
    "     - 1. Dataset and features description\n",
    "     - 2. Primary data analysis\n",
    "     - 3. Primary visual analysis of the features\n",
    "     - 4. Patterns, insights, pecularities of data\n",
    "     - 5. Metrics selection\n",
    "     - 6. Model selection\n",
    "     - 7. Data preprocessing\n",
    "     - 8. Cross-validation and adjustment of model hyperparameters\n",
    "     - 9. Creation of new features and description of this process\n",
    "     - 10. Plotting training and validation curves \n",
    "     - 11. Prediction for test or hold-out samples\n",
    "     - 12. Conclusions\n",
    "     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "KVrwRDk6FDsi"
   },
   "source": [
    ">>> ## Part 1. Dataset and features description"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "iC7u2raFrhYl"
   },
   "source": [
    "> ### **Abstract**: \n",
    "\n",
    "\n",
    "> Rossmann operates over 3,000 drug stores in 7 European countries. Currently, Rossmann store managers are tasked with predicting their daily sales for up to six weeks in advance. Store sales are influenced by many factors, including promotions, competition, school and state holidays, seasonality, and locality. With thousands of individual managers predicting sales based on their unique circumstances, the accuracy of results can be quite varied.\n",
    "\n",
    "> In their first Kaggle competition, Rossmann is challenging us to predict 6 weeks of daily sales for 1,115 stores located across Germany. Reliable sales forecasts enable store managers to create effective staff schedules that increase productivity and motivation.\n",
    " \n",
    "> https://www.kaggle.com/c/rossmann-store-sales\n",
    "\n",
    "\n",
    "We are provided with historical sales data for 1,115 Rossmann stores. The task is to forecast the \"Sales\" column for the test set. Note that some stores in the dataset were temporarily closed for refurbishment.\n",
    "\n",
    "> ### Files : \n",
    "\n",
    "Donwloaded from : https://www.kaggle.com/c/rossmann-store-sales/data\n",
    "\n",
    "train.csv - historical data including Sales\n",
    "\n",
    "test.csv - historical data excluding Sales\n",
    "\n",
    "sample_submission.csv - a sample submission file in the correct format\n",
    "\n",
    "store.csv - supplemental information about the stores\n",
    "\n",
    "> ### Some preliminary thoughts :\n",
    "\n",
    "A sales forecast is a tool that can help almost any company. Many companies rely on human forecasts that are not of a constant quality. Other companies use a standard tool that is not flexible enough to suit their needs. This competition gives researchers a possibility to create their own solution, that really improves business.\n",
    "\n",
    "Since data consists of sales based on date and time, we can build models with time series data, so it will serve us to extract date-time-features for further analysis. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ChaJ4BPsFDsj"
   },
   "source": [
    "> ###  Input variables :\n",
    "\n",
    "#### Data fields\n",
    "\n",
    "Self-explanatory fields :\n",
    "1. Date\n",
    "\n",
    "\n",
    "2. DayOfWeek\n",
    "       \n",
    "       \n",
    "The following fields have descriptions :\n",
    "\n",
    "3. Id - an Id that represents a (Store, Date) duple within the test set\n",
    "\n",
    "\n",
    "4. Store - a unique Id for each store\n",
    "\n",
    "\n",
    "5. Open - an indicator for whether the store was open: 0 = closed, 1 = open\n",
    "\n",
    "\n",
    "6. StateHoliday - indicates a state holiday. Normally all stores, with few exceptions, are closed on state holidays. Note that all schools are closed on public holidays and weekends. a = public holiday, b = Easter holiday, c = Christmas, 0 = None\n",
    "\n",
    "\n",
    "7. SchoolHoliday - indicates if the (Store, Date) was affected by the closure of public schools\n",
    "\n",
    "\n",
    "8. StoreType - differentiates between 4 different store models: a, b, c, d\n",
    "\n",
    "\n",
    "9. Assortment - describes an assortment level: a = basic, b = extra, c = extended\n",
    "\n",
    "\n",
    "10. CompetitionDistance - distance in meters to the nearest competitor store\n",
    "\n",
    "\n",
    "11. CompetitionOpenSince[Month/Year] - gives the approximate year and month of the time the nearest competitor was opened\n",
    "\n",
    "\n",
    "12. Promo - indicates whether a store is running a promo on that day\n",
    "\n",
    "\n",
    "13. Promo2 - Promo2 is a continuing and consecutive promotion for some stores: 0 = store is not participating, 1 = store is participating\n",
    "\n",
    "\n",
    "14. Promo2Since[Year/Week] - describes the year and calendar week when the store started participating in Promo2\n",
    "\n",
    "\n",
    "15. PromoInterval - describes the consecutive intervals Promo2 is started, naming the months the promotion is started anew. E.g. \"Feb,May,Aug,Nov\" means each round starts in February, May, August, November of any given year for that store\n",
    "\n",
    "\n",
    "> ### Output variable (desired target):\n",
    "16. Sales - the turnover for any given day \n",
    "\n",
    "\n",
    "We also have variable, that obviously high correlated with target:\n",
    "17. Customers - the number of customers on a given day\n",
    "\n",
    "> ### Evaluation:\n",
    "\n",
    "Submissions are evaluated on the Root Mean Square Percentage Error (RMSPE). \n",
    "<center>\n",
    "<img src=\"https://habrastorage.org/webt/kq/pz/hx/kqpzhxo3app_0yojqxkxa42uzx0.png\" />\n",
    "where ${y_i}$ denotes the sales of a single store on a single day and ${y(hat)_i}$ denotes the corresponding prediction. Any day and store with 0 sales is ignored in scoring.\n",
    "    \n",
    "Also it's worth saying: obviously there will be not any Sales in closed stores in particular days, so we can put 0 in the test target in rows, where store is closed. But since in the description of this competetion is said that Sales of closed stores will be ignored in resulting score - we won't pay attention to this rows in test data.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "_IFZI6k6FDsk"
   },
   "source": [
    ">>> ## Part 2. Primary data analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">### Loading packages :\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 17
    },
    "colab_type": "code",
    "id": "a2m9nQz3FDsm",
    "outputId": "9ce3aad7-1ad4-4b18-c870-8b4adee47481"
   },
   "outputs": [],
   "source": [
    "#system\n",
    "import os\n",
    "import sys\n",
    "import  gc\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "\n",
    "# basic\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from pandas import datetime\n",
    "import pickle\n",
    "import scipy\n",
    "import math\n",
    "import collections\n",
    "import sklearn\n",
    "\n",
    "\n",
    "# data preproc\n",
    "from sklearn.preprocessing import MinMaxScaler\n",
    "from scipy.stats import boxcox\n",
    "from scipy.special import inv_boxcox\n",
    "\n",
    "# data visualization\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "%matplotlib inline\n",
    "import itertools\n",
    "\n",
    "\n",
    "# statistics\n",
    "import scipy.stats as st\n",
    "from statsmodels.distributions.empirical_distribution import ECDF\n",
    "#from sklearn import metrics\n",
    "from sklearn.metrics import mean_squared_error, make_scorer\n",
    "\n",
    "\n",
    "# models\n",
    "import xgboost as xgb\n",
    "from sklearn.model_selection import GridSearchCV, train_test_split, RandomizedSearchCV, validation_curve\n",
    "from sklearn.model_selection import KFold, learning_curve\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.multioutput import MultiOutputClassifier\n",
    "from sklearn.ensemble import RandomForestClassifier\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### Importing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Path to our data\n",
    "PATH_TO_alldata = os.path.join(os.path.abspath(os.curdir), 'all_data')\n",
    "\n",
    "# importing train and test data\n",
    "initial_train = pd.read_csv(os.path.join(PATH_TO_alldata, \"train.csv\"), parse_dates = True, index_col = 'Date')\n",
    "train = initial_train.copy()\n",
    "initial_test = pd.read_csv(os.path.join(PATH_TO_alldata, \"test.csv\"), parse_dates = True, index_col = 'Date')\n",
    "test = initial_test.copy()\n",
    "\n",
    "# additional store data\n",
    "store = pd.read_csv(os.path.join(PATH_TO_alldata, \"store.csv\"))\n",
    "# As we did it with tran and test, create copy of inital frame and work with it:\n",
    "initial_store = store.copy()\n",
    "\n",
    "# reading sample_submission \n",
    "sample_submission = pd.read_csv(os.path.join(PATH_TO_alldata, \"sample_submission.csv\")) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### Main info about our data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"\\n----------------------------train ----------------------------\\n\")\n",
    "train.info()\n",
    "print(\"\\n----------------------------store ----------------------------\\n\")\n",
    "store.info()\n",
    "print(\"\\n----------------------------test ----------------------------\\n\")\n",
    "test.info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Number of Rows TRAIN:\", train.shape[0])\n",
    "print(\"Number of Columns TRAIN:\", train.shape[1])\n",
    "print(\"---------------------------\")\n",
    "print(\"Number of Rows TEST:\", test.shape[0])\n",
    "print(\"Number of Columns TEST:\", test.shape[1])\n",
    "print(\"---------------------------\")\n",
    "print(\"Number of Rows STORE:\", store.shape[0])\n",
    "print(\"Number of Columns STORE:\", store.shape[1])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### Let's look at heads :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "store.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some easy statistics: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### Counting number of NANs in each column of a Data Frame:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "def nans_counter(df):\n",
    "    for column in df.columns.values:\n",
    "        print(column,  \"NANs : \", pd.isnull(df.loc[:, column]).sum())\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"\\n TRAIN \")\n",
    "nans_counter(train)\n",
    "print(\"\\n TEST \")\n",
    "nans_counter(test)\n",
    "print(\"\\n STORE \")\n",
    "nans_counter(store)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we see that our train df doesn't have any NANs, test df has 11 NANs in OPEN column. \n",
    "And store df has many nans . \n",
    "We will consider  and remove this NANs later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### Now we will consider influence of different columns on target variable:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Firstly, we will take **Open** column and look on additional info about **SchoolHoliday**, **StateHoliday**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Over entire train period, {} - is the number of times that different stores closed on given days.\"\n",
    "       .format(train[(train.Open == 0)].count()[0]), \"\\n\")\n",
    "\n",
    "print(\"{} - is the number of closed stores because of a school holiday. \"\n",
    "       .format(train[(train.Open == 0) &\n",
    "                     (train.SchoolHoliday == 1) &\n",
    "                     (train.StateHoliday == '0')]\n",
    "               .count()[0]), \"\\n\")\n",
    "\n",
    "print(\"{} - is the number of closed stores  because of either a bank holiday or easter or christmas.\"\n",
    "       .format(train[(train.Open == 0) &\n",
    "                     ((train.StateHoliday == 'a') |\n",
    "                      (train.StateHoliday == 'b') | \n",
    "                      (train.StateHoliday == 'c'))]\n",
    "               .count()[0]), \"\\n\")\n",
    "\n",
    "print(\"{} - is the number of times, when shops were closed on days for no apparent reason when no holiday was announced.\"\n",
    "      .format(train[(train.Open == 0) &\n",
    "                    (train.StateHoliday == \"0\") &\n",
    "                    (train.SchoolHoliday == 0)]\n",
    "              .count()[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the descrition of the this task we can get, that Rossman stores had undergoing refurbishments sometimes and has been closed. Most probably those times, when stores were closed without any reason (in dataset), mean this refurbishments. \n",
    "\n",
    "And since we don't want to bias our decision, wi will not consider those exceptions, so we will get rid of closed stores and prevent the models to train on them and get false guidance.\n",
    "\n",
    "In this case we will analyse only open stores since a close store yield a profit of 0.\n",
    "And as I said in the beginning: we will not pay any attention to similar situations in test data. \n",
    "\n",
    "And distribution of **Opened-closed stores per day of week** is here : \n",
    "> ### 'Open' option with respect to Day Of week:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axis = plt.subplots(1,1,figsize=(15,4))\n",
    "sns.countplot(x='Open',hue='DayOfWeek', data=train,\n",
    "              palette=sns.color_palette(\"Set1\", n_colors=7), ax=axis)\n",
    "plt.title(\"Train Data\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# getting rid of closed stores\n",
    "train = train[(train.Open != 0) & (train.Sales != 0)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we remember: we have **11 NANs in Open column in Test Data**. \n",
    "\n",
    "So it's time to deal with it. First of all, we have to know, which store has this NANs. Then we will explore: when this store was closed. After that we will decide whether every NAN should have a replacing option = 0 or not. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axis = plt.subplots(1,1,figsize=(15,4))\n",
    "sns.countplot(x='Open',hue='DayOfWeek', data=test,\n",
    "              palette=sns.color_palette(\"Set1\", n_colors=7), ax=axis)\n",
    "plt.title(\"Test Data\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test[(test['Open'] != 0) & (test['Open']!=1)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Ok, that's only 622 store, that has NANs. Let's consider days with NAN :\\n\", \n",
    "      sorted(test[(test['Open'] != 0) & (test['Open']!=1)].DayOfWeek.unique()))\n",
    "\n",
    "fig, axis = plt.subplots(1,1,figsize=(15,4))\n",
    "sns.countplot(x='Open',hue='DayOfWeek', data=test[test['Store']==622],\n",
    "              palette=sns.color_palette(\"Set1\", n_colors=7), ax=axis)\n",
    "plt.title(\"Test Data for store 622\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Paying attention to plot with test data and open column for 622 store-> we see, that on days [1..6] this store always open. \n",
    "\n",
    "It leads us to **put 1 in all NANs in this data**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# replacing NANs with OPEN-state\n",
    "test.Open.fillna(1, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now our data without useless information about closed stores. \n",
    "\n",
    "We can group by **date** and get **average sales**, and **percentage change** between the current and a prior element:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "average_sales    = train.groupby('Date')[\"Sales\"].mean()\n",
    "pct_change_sales = train.groupby('Date')[\"Sales\"].sum().pct_change()\n",
    "\n",
    "fig, (axis1,axis2) = plt.subplots(2,1, sharex=True, figsize=(15,8))\n",
    "\n",
    "# average sales\n",
    "ax1 = average_sales.plot(legend=True, ax=axis1, marker='p', title=\"Average Sales\")\n",
    "ax1.set_xticks(range(len(average_sales)))\n",
    "ax1.set_xticklabels(average_sales.index.tolist(), rotation=90)\n",
    "\n",
    "# percent change for sales\n",
    "ax2 = pct_change_sales.plot(legend=True,ax=axis2, marker='o', rot=90, colormap=\"summer\", title=\"Sales Percent Change\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Got from this plots: clearly seeing increasing of Sales in period around NewYear, seems to be real. \n",
    "\n",
    "Pattern of every year Average Sales Function is repeated every year - Ok, nothing strange, nothing super interesting. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">  ### So what about distribution of  **Sales** in different **days of week**?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(131)\n",
    "sns.barplot(x=train.groupby('DayOfWeek')['Sales'].sum().index, y=train.groupby('DayOfWeek')['Sales'].sum().values)\n",
    "ax1.set_title(\"Summing sales over each day\")\n",
    "\n",
    "ax2 = fig.add_subplot(132)\n",
    "sns.barplot(x=train.groupby('DayOfWeek')['Sales'].mean().index, y=train.groupby('DayOfWeek')['Sales'].mean().values)\n",
    "ax2.set_title(\"Mean sales in each day\")\n",
    "\n",
    "ax3 = fig.add_subplot(133)\n",
    "sns.countplot(train.DayOfWeek)\n",
    "ax3.set_title(\"Total count of sales in each day\")\n",
    "\n",
    "plt.rcParams['figure.figsize'] = 18, 7\n",
    "sns.set(font_scale=1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.boxplot(x='DayOfWeek', y='Sales', data=train, palette=\"muted\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, from plots we can see that highest profit is happened at 1th day.\n",
    "\n",
    "At 7th day there is very small quantity of sales (most probably because of closed stores), but mean-sale in this day one of the highest! There are 2 leaders in \"mean measure\" : 1th and 7th days.\n",
    "\n",
    "Also we can consider such interesting variable as **Promo** (whether a store is running a promo on that day or not). \n",
    "\n",
    " \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### Here we can see how depends **sales** and **customers** on every **day of week** and whether **promo** is on day or not :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"\\t\\t\\t----------Dependency of Sales on DayOfWeek with/without Promo--------\\n\")\n",
    "print(\"1nd bar is for day with Promo, 2nd bar is for day without Promo\")\n",
    "fig = plt.figure()\n",
    "# Bar width\n",
    "width = 0.4\n",
    "\n",
    "ax1 = fig.add_subplot(121)\n",
    "Day_Promo = train[train['Promo'] == 1].groupby('DayOfWeek')['Sales'].sum()\n",
    "Day_NoPromo = train[train['Promo'] != 1].groupby('DayOfWeek')['Sales'].sum()\n",
    "Day_Promo.plot(kind='bar', width=width, color='#FF0000', position=1, label='Day_Promo', rot=0)\n",
    "Day_NoPromo.plot(kind='bar', width=width, color='#0000FF', position=0, label='Day_NoPromo', rot=0)\n",
    "plt.xlabel('Day of week')\n",
    "plt.ylabel('sum')\n",
    "ax1.set_title(\"Summing Sales over each day\")\n",
    "\n",
    "ax2 = fig.add_subplot(122)\n",
    "Day_Promo = train[train['Promo'] == 1].groupby('DayOfWeek')['Sales'].mean()\n",
    "Day_NoPromo = train[train['Promo'] != 1].groupby('DayOfWeek')['Sales'].mean()\n",
    "Day_Promo.plot(kind='bar', width=width, color='#FF0000', position=1, label='Day_Promo', rot=0)\n",
    "Day_NoPromo.plot(kind='bar', width=width, color='#0000FF', position=0, label='Day_NoPromo', rot=0)\n",
    "plt.xlabel('Day of week')\n",
    "plt.ylabel('mean')\n",
    "ax2.set_title(\"Mean Sales in each day\")\n",
    "\n",
    "plt.rcParams['figure.figsize'] = 18, 8\n",
    "sns.set(font_scale=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"\\t\\t\\t----------Dependency of Customers on DayOfWeek with/without Promo--------\\n\")\n",
    "print(\"1nd bar is for day with Promo, 2nd bar is for day without Promo\")\n",
    "fig = plt.figure()\n",
    "# Bar width\n",
    "width = 0.4\n",
    "\n",
    "ax1 = fig.add_subplot(121)\n",
    "Day_Promo = train[train['Promo'] == 1].groupby('DayOfWeek')['Customers'].sum()\n",
    "Day_NoPromo = train[train['Promo'] != 1].groupby('DayOfWeek')['Customers'].sum()\n",
    "Day_Promo.plot(kind='bar', width=width, color='#FF0000', position=1, label='Day_Promo', rot=0)\n",
    "Day_NoPromo.plot(kind='bar', width=width, color='#0000FF', position=0, label='Day_NoPromo', rot=0)\n",
    "plt.xlabel('Day of week')\n",
    "plt.ylabel('sum')\n",
    "ax1.set_title(\"Summing Customers over each day\")\n",
    "\n",
    "ax2 = fig.add_subplot(122)\n",
    "Day_Promo = train[train['Promo'] == 1].groupby('DayOfWeek')['Customers'].mean()\n",
    "Day_NoPromo = train[train['Promo'] != 1].groupby('DayOfWeek')['Customers'].mean()\n",
    "Day_Promo.plot(kind='bar', width=width, color='#FF0000', position=1, label='Day_Promo', rot=0)\n",
    "Day_NoPromo.plot(kind='bar', width=width, color='#0000FF', position=0, label='Day_NoPromo', rot=0)\n",
    "plt.xlabel('Day of week')\n",
    "plt.ylabel('mean')\n",
    "ax2.set_title(\"Mean Customers in each day\")\n",
    "\n",
    "plt.rcParams['figure.figsize'] = 18, 8\n",
    "sns.set(font_scale=1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Got some interesting results here. Sales on saturday without Promo really tends to peak a maximum. \n",
    "On sunday, I think, there are many closed stores - so minimum benefits. Other days : we can see, that sales on Promo-days significantlly bigger.\n",
    "\n",
    "Interesting, that \"mean-measure\" of sales in Promo-Sunday doesn not have big difference from other days, but the total count of customers of this day is extralarge, and it leads total sum of sales to be bigger. \n",
    "\n",
    "In the beginning I was talking about obvious assumption about strong connection between **Sales** and **Customers**, now it's time to prove it precisely :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Pearson correlation : \", round(st.pearsonr(train.Customers, train.Sales)[0], 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our pearson correlation factor of 0.824 explains that there is a strong positive correlation between Sales and Customers. In general, the more customers you have in a store, the higher your sales for the day.\n",
    " \n",
    "\n",
    "So it seems interesting to create new variable : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# combinig new variable\n",
    "train['SalePerCustomer'] = train['Sales'] / train['Customers']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### Distribution analysis of target variables:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creating some usefull functions for distribution analysis: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def norm_testing(vector):\n",
    "    print(\"Skewness value :\", np.round(st.skew(vector.values),decimals=2))\n",
    "    print(\"The kurtosis for Fisher’s definition  :\", np.round(st.kurtosis(vector.values) , decimals=2))\n",
    "    print(\"Assymetry: \", np.abs(np.round(st.skew(vector.values),decimals=2)))\n",
    "    print(\"Excess: \", np.abs(np.round(st.kurtosis(vector.values) ,decimals=2)))\n",
    "    \n",
    "    print(\"The kurtosis test that it's a normal distribution :\", \n",
    "         \"\\nZ-score: {}  \\nThe 2-sided p-value for the hypothesis test: {}\".\n",
    "         format(round(st.kurtosistest(vector.values)[0], 4), st.kurtosistest(vector.values)[1]))\n",
    "    \n",
    "    # for getting better accuracy in Sapiro-Wilk test, let's take from data, for example, 4500 random values\n",
    "    # In documentation we can find max quantity of values for accurate results : 5000 \n",
    "    vector = vector.reset_index(drop=True) #reindexing to [0,1,2,3 etc.] \n",
    "    indices = np.random.choice(a=vector.index, size=4500, replace=False)\n",
    "    \n",
    "    print(\"The Shapiro-Wilk test for the null hypothesis, that it's a normal distribution: \",\n",
    "          \"\\nThe test statistic: {} \\nThe p-value for the hypothesis test {}\".\n",
    "          format(round(st.shapiro(vector[indices].values)[0], 4), st.shapiro(vector[indices].values)[1]))\n",
    "    \n",
    "\n",
    "def get_distrib_stats(data):\n",
    "    norm_testing(data)\n",
    "    \n",
    "    sns.set(rc={\"figure.figsize\": (20, 5)})\n",
    "    sns.set(font_scale=1.2)\n",
    "    fig = plt.figure()\n",
    "    ax1 = fig.add_subplot(121)\n",
    "    p_x = sns.distplot(data, fit=st.norm, kde=True,ax=ax1, bins=50)\n",
    "    #ax1.legend()\n",
    "    ax1.set_title(\"Probability density function\")\n",
    "    ax2 = fig.add_subplot(122)\n",
    "    prob = st.probplot(data, dist=st.norm, plot=ax2)\n",
    "    ax2.set_title('Probabilyty plot')\n",
    "    #ax2.legend()\n",
    "    plt.show()\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"\\t\\t\\t-----------Distribution analysis for Sales-----------\")\n",
    "get_distrib_stats(train['Sales'])\n",
    "\n",
    "print(\"\\t\\t\\t-----------Distribution analysis for Customers-----------\")\n",
    "get_distrib_stats(train['Customers'])\n",
    "\n",
    "print(\"\\t\\t\\t-----------Distribution analysis for SalePerCustomer-----------\")\n",
    "get_distrib_stats(train['SalePerCustomer'])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### OK, it's not normally distributed. \n",
    "\n",
    "Better for our future modeling to have target variable normally distributed.\n",
    "\n",
    "How we can deal with it : we can remove \"unconvienient\" values, like large tails, so it gives us a bit different picture: data seems to become closer to normal. But here we have to deal with some data leaking.\n",
    "Also we can add log() function to our target and we will gain again some \"improvment\" in the way getting variable more normally distributed. \n",
    "\n",
    "For sure we can consider more different techniques, I tried many, and for \"place and time\" saving I will show the best way I explored this data to make variable normally distributed.\n",
    "\n",
    "We will use **coxbox**.\n",
    "\n",
    "But first let's look at **boxplots** of Sales and Customers: \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(121)\n",
    "sns.violinplot(x='Sales', scale='count',data=train[train['Sales'] > 0])\n",
    "ax1.set_title(\"violinplot of sales\")\n",
    "\n",
    "ax2 = fig.add_subplot(122)\n",
    "sns.boxplot( x='Sales', data=train[train['Sales'] > 0], orient=\"h\")\n",
    "ax2.set_title(\"Boxplot of sales\")\n",
    "\n",
    "plt.rcParams['figure.figsize'] = 18, 7\n",
    "sns.set(font_scale=3)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(121)\n",
    "sns.violinplot(x='Customers', scale='count',data=train[train['Customers'] > 0])\n",
    "ax1.set_title(\"violinplot of Customers\")\n",
    "\n",
    "ax2 = fig.add_subplot(122)\n",
    "sns.boxplot( x='Customers', data=train[train['Customers'] > 0], orient=\"h\")\n",
    "ax2.set_title(\"Boxplot of Customers\")\n",
    "\n",
    "plt.rcParams['figure.figsize'] = 18, 7\n",
    "sns.set(font_scale=1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are some outliers in target variables, so we will decide how to deal with it in next chapters. \n",
    "\n",
    "And now **coxboxing the target** : \n",
    "\n",
    "(we will also save **lambdas** from coxbox function for next usage in **inverse coxbox** function later, when prediction will be calculated, and so we will get true prediction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# coxboxing and save lambdas : \n",
    "lambdas_coxbox = []\n",
    "\n",
    "coxboxed, lam = boxcox(train['Sales'].values)\n",
    "sales_coxboxed = pd.Series(coxboxed, index=train.Sales.index)\n",
    "lambdas_coxbox.append(lam)\n",
    "\n",
    "coxboxed, lam = boxcox(train['Customers'].values)\n",
    "cust_coxboxed = pd.Series(coxboxed, index=train.Customers.index)\n",
    "lambdas_coxbox.append(lam)\n",
    "\n",
    "coxboxed, lam = boxcox(train['SalePerCustomer'].values)\n",
    "sale_per_cust_coxboxed = pd.Series(coxboxed, index=train.SalePerCustomer.index)\n",
    "lambdas_coxbox.append(lam)\n",
    "\n",
    "print('Lambda for Sales: %f' % lambdas_coxbox[0])\n",
    "print('Lambda for Customers: %f' % lambdas_coxbox[1])\n",
    "print('Lambda for SalePerCustomer: %f' % lambdas_coxbox[2])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"\\t\\t\\t-----------Distribution analysis for coxboxed Sales-----------\")\n",
    "get_distrib_stats(sales_coxboxed)\n",
    "\n",
    "print(\"\\t\\t\\t-----------Distribution analysis for coxboxed Customers-----------\")\n",
    "get_distrib_stats(cust_coxboxed)\n",
    "\n",
    "print(\"\\t\\t\\t-----------Distribution analysis for coxboxed SalePerCustomer-----------\")\n",
    "get_distrib_stats(sale_per_cust_coxboxed)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, now it seems to be better. \n",
    "Also, we didn't do anything with outliers, so, perhaps, removing this values will lead us to better distribuation without leaking in accuracy of forecasting.\n",
    "\n",
    "Let's pay attention last time for now on distribution statistics and also check, how much data we will have after **removing outliers** at all? : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# finding outliers\n",
    "fig = plt.figure()\n",
    "\n",
    "ax1 = fig.add_subplot(131)\n",
    "sns.boxplot(data=sales_coxboxed, orient=\"h\")\n",
    "ax1.set_title(\"Boxplot of coxboxed sales\")\n",
    "\n",
    "ax2 = fig.add_subplot(132)\n",
    "sns.boxplot(data=cust_coxboxed, orient=\"h\")\n",
    "ax2.set_title(\"Boxplot of coxboxed Customers\")\n",
    "\n",
    "ax3 = fig.add_subplot(133)\n",
    "sns.boxplot(data=sale_per_cust_coxboxed, orient=\"h\")\n",
    "ax3.set_title(\"Boxplot of coxboxed SalePerCustomer\")\n",
    "\n",
    "\n",
    "\n",
    "plt.rcParams['figure.figsize'] = 17, 14\n",
    "sns.set(font_scale=1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"\\t\\t\\t-----------Distribution analysis for cutted coxboxed Sales-----------\")\n",
    "get_distrib_stats(sales_coxboxed.where((sales_coxboxed > 8) & (sales_coxboxed < 15.5 ))\n",
    "                            .sort_values().dropna())\n",
    "\n",
    "print(\"\\t\\t\\t-----------Distribution analysis for cutted coxboxed Customers-----------\")\n",
    "get_distrib_stats(cust_coxboxed.where((cust_coxboxed > 3.1) & (cust_coxboxed < 4.5))\n",
    "                            .sort_values().dropna())\n",
    "\n",
    "print(\"\\t\\t\\t-----------Distribution analysis for cutted coxboxed SalePerCustomer-----------\")\n",
    "get_distrib_stats(sale_per_cust_coxboxed.where(sale_per_cust_coxboxed < 6)\n",
    "                            .sort_values().dropna())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"We will save :\")\n",
    "print(\"{0:.3f}% of data for Sales\".format(len(sales_coxboxed.where((sales_coxboxed > 8) & (sales_coxboxed < 15.5 ))\n",
    "                            .sort_values().dropna()) / len(sales_coxboxed.values)*100))\n",
    "\n",
    "print(\"{0:.3f}% of data for Customers\".format(len(cust_coxboxed.where((cust_coxboxed > 3.1) & (cust_coxboxed < 4.5))\n",
    "                            .sort_values().dropna()) / len(cust_coxboxed.values)*100))\n",
    "\n",
    "print(\"{0:.3f}% of data for SalePerCustomer\".format(len(sale_per_cust_coxboxed.where(sale_per_cust_coxboxed < 6)\n",
    "                            .sort_values().dropna()) / len(sale_per_cust_coxboxed.values)*100))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we get pretty nice results, we are going to change values for all three variables in train data set, and then we will go to next part of out research: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# take coxboxed values\n",
    "train.Sales = sales_coxboxed\n",
    "train.Customers = cust_coxboxed\n",
    "train.SalePerCustomer = sale_per_cust_coxboxed\n",
    "\n",
    "# removing outliers from each variable\n",
    "print(\"rows before cutting : \", train.shape[0])\n",
    "train = train[(train.Sales > 8) & (train.Sales < 15.5 )]\n",
    "train = train[(train.Customers > 3.1) & (train.Customers < 4.5 )]\n",
    "train = train[train.SalePerCustomer < 6]\n",
    "print(\"rows after cutting : \", train.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">>> ## Part 3 & Part 4. Primary visual data analysis &  Patterns, insights, peculiarities of data "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have spent a lot of time on looking for interesting moments in target variables, and now we will start from creating new self-explanatory features and showing visual analysis.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train['Year'] = train.index.year\n",
    "train['Month'] = train.index.month\n",
    "train['Day'] = train.index.day\n",
    "train['WeekOfYear'] = train.index.weekofyear\n",
    "\n",
    "test['Year'] = test.index.year\n",
    "test['Month'] = test.index.month\n",
    "test['Day'] = test.index.day\n",
    "test['WeekOfYear'] = test.index.weekofyear\n",
    "\n",
    "#print(\"\\n TRAIN \")\n",
    "#nans_counter(train)\n",
    "#print(\"\\n TEST \")\n",
    "#nans_counter(test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": false
   },
   "source": [
    "Since we started, we didn't dive deep inside into store dataframe. \n",
    "\n",
    "So it's time do it: \n",
    "\n",
    "> ### Diving deep inside into store dataframe: \n",
    "\n",
    "Firstly, as I said before, we will deal with missing values. \n",
    "\n",
    "Then we will merge our train and test frames with information in store-df. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We have already explored about nans : \n",
    " **STORE**\n",
    " \n",
    "Store NANs :  0\n",
    "\n",
    "StoreType NANs :  0\n",
    "\n",
    "Assortment NANs :  0\n",
    "\n",
    "CompetitionDistance NANs :  3\n",
    "\n",
    "CompetitionOpenSinceMonth NANs :  354\n",
    "\n",
    "CompetitionOpenSinceYear NANs :  354\n",
    "\n",
    "Promo2 NANs :  0\n",
    "\n",
    "Promo2SinceWeek NANs :  544\n",
    "\n",
    "Promo2SinceYear NANs :  544\n",
    "\n",
    "PromoInterval NANs :  544\n",
    "\n",
    "\n",
    "Ok, pretty obvious, that **Promo2SinceWeek**, **Promo2SinceYear** and **PromoInterval** have NANs because of zero values in Promo2: in meaning of absent continuous promotional activities for those stores. Since we understand this fact, we will just replace NANs in this columns with 0 : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "store.Promo2SinceWeek.fillna(0,inplace=True)\n",
    "store.Promo2SinceYear.fillna(0,inplace=True)\n",
    "store.PromoInterval.fillna(0,inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Second problem - is **CompetitionDistance** NANs.\n",
    "\n",
    "Let's explore it values and distribution : \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "store_1 = store.drop(store[pd.isnull(store.CompetitionDistance)].index)\n",
    "fig, axes = plt.subplots(1, 2, figsize=(17,3.5))\n",
    "\n",
    "axes[0].boxplot(store_1.CompetitionDistance, showmeans=True, vert=False)\n",
    "axes[0].set_title('Boxplot for Closest Competition')\n",
    "\n",
    "axes[1].hist(store_1.CompetitionDistance, bins=50)\n",
    "axes[1].set_title(\"Closest Competition histogram\")\n",
    "\n",
    "# deleting redundant Data Frame\n",
    "del store_1\n",
    "gc.collect()\n",
    "\n",
    "print(\"Mean : \", np.nanmean(store.CompetitionDistance), \n",
    "      \"\\nMedian : \", np.nanmedian(store.CompetitionDistance), \n",
    "      \"\\nStandard Dev : \", np.nanstd(store.CompetitionDistance))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a highly skewed distribution for this variable with a huge difference between the mean and the median. This being caused by the significant value of dispersion with a standard deviation of 7659, higher than the mean and the median. \n",
    "\n",
    "So unlike mean, the median gives us the positional center of the data so it is not influenced by outliers. And we will replace NANs with the median.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "store['CompetitionDistance'].fillna(store['CompetitionDistance'].median(), inplace = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last problem here - is NANs in **CompetitionOpenSinceMonth** and **CompetitionOpenSinceYear**.\n",
    "\n",
    "The easiest way here is fullfil the fields with 0. \n",
    "\n",
    "More interesting way: take top Month and Year from distributions of this values in dataset, but it can lead to wrong results because of this \"inaccurate\" transformation. \n",
    "\n",
    "And the most creative way could be to apply a **multilabel classification** algorithm and train on the non Nan fields and then predict what could be most probably the month and year for those fields. \n",
    "(this idea I explored while reading the discussion and forums for this competition. I didn't find in kernels, that anyone made something like this, so this mini-task is going to be more intriguing)\n",
    "\n",
    "> So we will try to do it with very basic algorithm and super easy feature auto selection,  because of the main purpose of this work: we are going to do research in another forecasting task, not multilabel classification. But it's interesting to do fast implementation for such NANs replacing.  \n",
    "\n",
    "I think, it's better to do full implementation of this idea in function, returning non NAN columns, because of memory saving. \n",
    "\n",
    "> ### Function for multilabel classification prediction inplace NANs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mlc = multilabel classification\n",
    "def mlc_fun(df, trees=1000):\n",
    "    # ##### split on train and predict sets\n",
    "    # taking only not NAN objects in train : \n",
    "    df_store_tr = df[df['CompetitionOpenSinceMonth'] > 0].copy()\n",
    "    \n",
    "    # dropping all not NAN objects from df for test:     \n",
    "    df_store_te = df.copy()\n",
    "    df_store_te.drop(df_store_tr.index, inplace=True, axis=0)\n",
    "    print(\">>>checking correction of separating on train and test: \\n\",\n",
    "          df_store_tr.shape[0], '+', df_store_te.shape[0], \"=\",\n",
    "          df_store_tr.shape[0] + df_store_te.shape[0], \"== full length:\", df.shape[0])\n",
    "        \n",
    "    # taking store values for merge and the end\n",
    "    store_tr = df_store_tr.Store.values\n",
    "    store_te = df_store_te.Store.values\n",
    "    \n",
    "    # taking true target varibales for train. And store column for merge dataframes later \n",
    "    y_1_2_true = df_store_tr[['CompetitionOpenSinceMonth', 'CompetitionOpenSinceYear', 'Store']]\n",
    "    \n",
    "    # we will not use Store column in prediction - so drop it. Also drop target columns\n",
    "    df_store_tr.drop(columns=['CompetitionOpenSinceMonth', 'CompetitionOpenSinceYear', 'Store'], inplace=True)\n",
    "    df_store_te.drop(columns=['CompetitionOpenSinceMonth', 'CompetitionOpenSinceYear', 'Store'], inplace=True)\n",
    "    \n",
    "    # ##### separating target for mini-preprocessing\n",
    "    y_1 = y_1_2_true.CompetitionOpenSinceMonth\n",
    "    y_2 = y_1_2_true.CompetitionOpenSinceYear\n",
    "    \n",
    "    # taking unique monthes that exist in train set\n",
    "    monthes_ = [\"month{}\".format(i) for i in sorted(y_1.unique())]\n",
    "    # dict for replace values in target with \"str\"type ones\n",
    "    d1 = {k: v for (k,v) in zip(range(1,13), monthes_)}\n",
    "    # inverse dict for de-code the predicted target in future \n",
    "    d1_ = {v: k for (k,v) in zip(range(1,13), monthes_)}  \n",
    "    y_1 = y_1.map(d1)\n",
    "    \n",
    "    # taking unique years that exist in train set\n",
    "    years_ = [\"year{}\".format(i) for i in sorted(y_2.unique())]\n",
    "    # dict for replace values in target with \"str\"type ones\n",
    "    d2 = {k: v for (k,v) in zip(sorted(y_2.unique()), years_)}\n",
    "    # inverse dict for de-code the predicted target in future \n",
    "    d2_ = {v: k for (k,v) in zip(sorted(y_2.unique()), years_)}\n",
    "    y_2 = y_2.map(d2)\n",
    "\n",
    "    # full stack target vectors \n",
    "    Y_ = np.vstack((y_1, y_2)).T\n",
    "    # memory cleaning \n",
    "    del y_1, y_2\n",
    "    gc.collect()\n",
    "\n",
    "    # ##### super fast feature generating and coding\n",
    "    print(\">>>getting dummies\")\n",
    "    # union of train and test for getting dummies\n",
    "    train_objs_num = len(df_store_tr)\n",
    "    dataset = pd.concat(objs=[df_store_tr, df_store_te], axis=0)\n",
    "    dataset_preprocessed = pd.get_dummies(dataset)\n",
    "    df_store_tr_dum = dataset_preprocessed[:train_objs_num]\n",
    "    df_store_te_dum = dataset_preprocessed[train_objs_num:]\n",
    "    # memory cleaning \n",
    "    del dataset_preprocessed, dataset, train_objs_num\n",
    "    gc.collect()\n",
    "    \n",
    "    #df_store_tr_dum = pd.get_dummies(df_store_tr)\n",
    "    #df_store_te_dum = pd.get_dummies(df_store_te)\n",
    "    \n",
    "    # ##### implementing easy model for classification ot multi target\n",
    "    print(\">>>starting classification model\")\n",
    "    forest = RandomForestClassifier(n_estimators=trees, random_state=17)\n",
    "    multi_target_forest = MultiOutputClassifier(forest, n_jobs=4)\n",
    "    y_pred = multi_target_forest.fit(df_store_tr_dum, Y_).predict(df_store_te_dum)\n",
    "    \n",
    "    # ##### forming convinient view of predicted values\n",
    "    print(\">>>prediction Done. \\nReturning...\")\n",
    "    y_pred_ = df_store_te_dum.copy()\n",
    "    y_pred_['CompetitionOpenSinceMonth'] = y_pred.T[0]\n",
    "    y_pred_['CompetitionOpenSinceYear'] = y_pred.T[1]\n",
    "    y_pred_.drop(columns=df_store_te_dum.columns, inplace=True)\n",
    "    # getting back store column\n",
    "    y_pred_['Store'] = store_te\n",
    "    \n",
    "    # decoding with inverse dicts\n",
    "    y_pred_.CompetitionOpenSinceMonth = y_pred_.CompetitionOpenSinceMonth.map(d1_)\n",
    "    y_pred_.CompetitionOpenSinceYear = y_pred_.CompetitionOpenSinceYear.map(d2_)\n",
    "    \n",
    "    # concatinate true target and predicted\n",
    "    full_target = pd.concat([y_pred_, y_1_2_true])\n",
    "    \n",
    "    # ##### creating new dataframe for returning option\n",
    "    df_store_new_1 = df.copy()\n",
    "    df_store_new_1.drop(columns=['CompetitionOpenSinceMonth', 'CompetitionOpenSinceYear'], inplace=True)\n",
    "    df_store_new_2 = pd.merge(df_store_new_1, full_target, how = 'left', on = 'Store')\n",
    "    \n",
    "    # memory cleaning \n",
    "    del df_store_new_1, full_target, y_pred_, y_1_2_true, d1 ,d2, d1_, d2_\n",
    "    del df_store_tr, df_store_te, store_tr, store_te, df_store_tr_dum, df_store_te_dum\n",
    "    del forest, multi_target_forest\n",
    "    gc.collect()\n",
    "    \n",
    "    return(df_store_new_2)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, function is ready. Let's just check it's accuracy of this prediction function on train data. \n",
    "\n",
    "For sure, it won't be high quality prediction, but I think, that if it is better than 0% - the initial idea of implementing classification wasn't for nothing! Because if we put just zeros inplace of NANs - it's around zero significance.\n",
    "\n",
    "So the statement is our **multilabel classification function has some importance!**\n",
    "\n",
    "For proving this this we will take existing values and replace them with 0, then we will compare results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# take only existing values\n",
    "store_check = store[store['CompetitionOpenSinceMonth'] > 0].copy()\n",
    "\n",
    "# random choicing values that will be put 0\n",
    "indc = store_check.index\n",
    "indc_r = np.random.choice(a=indc, size=int(len(indc)*0.2), replace=False)\n",
    "\n",
    "# saving true target\n",
    "y_true = store_check.loc[indc_r, ['CompetitionOpenSinceMonth', 'CompetitionOpenSinceYear']]\n",
    "\n",
    "# creating \"pseudo NANs\"\n",
    "store_check.loc[indc_r, 'CompetitionOpenSinceMonth'] = 0 \n",
    "store_check.loc[indc_r, 'CompetitionOpenSinceYear'] = 0\n",
    "\n",
    "# get prediction\n",
    "store_check_new = mlc_fun(store_check, 10000)\n",
    "\n",
    "# checking for evaluation\n",
    "print(\"for CompetitionOpenSinceMonth : \")\n",
    "check = y_true.CompetitionOpenSinceMonth == store_check_new.loc[indc_r, 'CompetitionOpenSinceMonth']\n",
    "check_count = check.apply(lambda x: 1 if x==True else 0)\n",
    "print(\"rigth\", sum(check_count), \"of\", len(indc_r), \"=> accuracy =\", round(sum(check_count)/len(indc_r)*100, 2), \"%\")\n",
    "\n",
    "print(\"for CompetitionOpenSinceYear: \")\n",
    "check = y_true.CompetitionOpenSinceYear == store_check_new.loc[indc_r, 'CompetitionOpenSinceYear']\n",
    "check_count = check.apply(lambda x: 1 if x==True else 0)\n",
    "print(\"rigth\", sum(check_count), \"of\", len(indc_r), \"=> accuracy =\", round(sum(check_count)/len(indc_r)*100, 2), \"%\")\n",
    "\n",
    "# cleaning memory\n",
    "del store_check, store_check_new, indc, indc_r, check, check_count\n",
    "gc.collect()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, we got **greater then 0 accuracy prediction!** Alright, it's super small for good classification, but we will go with this implemention and predict \"classes\" for true NANs.\n",
    "\n",
    "I think, this idea is really good, and here we have large space for make something creative and wellpredicting. But for this step, it's enough. We can always back here and change something for model improvment.\n",
    "\n",
    "And now we are going to the next point.\n",
    "\n",
    "It would better not rewrite the store-df. Perhaps, we will use both in some model later and check which one better. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# random forest with 10k estimators \n",
    "store_mlc = mlc_fun(store, 10000)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save it if you want\n",
    "# store_mlc.to_csv(os.path.join(PATH_TO_alldata, \"store_mlc.csv\"), index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# checking is there any NAN\n",
    "nans_counter(store_mlc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now it's clean, no NANs in store-df! \n",
    "\n",
    "Since train and store DFs are ready for every kind of analysis, we will start from merging them!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train_store_mlc = pd.merge(train, store_mlc, how = 'left', on = 'Store')\n",
    "\n",
    "# all the time we make similar actions with test data\n",
    "df_test_store_mlc = pd.merge(test, store_mlc, how = 'left', on = 'Store')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### Analysis Sales with respect to CompetitionDistance :\n",
    "\n",
    "I think, it's important charactristic, supposing, that stores without near competitors have significantly more sales then other stores.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis_comp_dist = df_train_store_mlc.loc[:, ['CompetitionDistance', 'Sales', 'Customers']].copy()\n",
    "analysis_comp_dist.CompetitionDistance.hist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, let's split CompetitionDistance for 5-6 parts and then look on AVERAGE of Sales and Customers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analysis_comp_dist['CompetitionDist_cut'] = pd.cut(analysis_comp_dist['CompetitionDistance'], 5)\n",
    "print(\"for Sales: \", analysis_comp_dist.groupby(by=\"CompetitionDist_cut\").Sales.mean())\n",
    "print(\"Customers: \", analysis_comp_dist.groupby(by=\"CompetitionDist_cut\").Customers.mean())\n",
    "\n",
    "del analysis_comp_dist\n",
    "gc.collect()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can state that the stores that are the furthest have the highest average sales and number of customers. \n",
    "\n",
    "But it is not so significant difference from the stores that have their competitors not so far.\n",
    "\n",
    "Anyway, if there are no competition nearby, stores tend to sell more and have more customers because there are almost a monopoly in this region. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### Analysis Sales, Customers and SalePerCustomer with respect to Store Type :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 3,figsize=(17,10))\n",
    "palette = itertools.cycle(sns.color_palette(\"Set1\", n_colors=4))\n",
    "clrs = [next(palette),next(palette),next(palette),next(palette)]\n",
    "\n",
    "axes[0,0].bar(store_mlc.groupby(by=\"StoreType\").count().Store.index, \n",
    "              store_mlc.groupby(by=\"StoreType\").count().Store, color=clrs)\n",
    "axes[0,0].set_title(\"Number of Stores per Store Type \\n Fig 1.1\")\n",
    "\n",
    "axes[0,1].bar(df_train_store_mlc.groupby(by=\"StoreType\").sum().Sales.index,\n",
    "              df_train_store_mlc.groupby(by=\"StoreType\").sum().Sales/1e9, color=clrs)\n",
    "axes[0,1].set_title(\"Total Sales per Store Type (in Billions) \\n Fig 1.2\")\n",
    "\n",
    "axes[0,2].bar(df_train_store_mlc.groupby(by=\"StoreType\").sum().Customers.index,\n",
    "              df_train_store_mlc.groupby(by=\"StoreType\").sum().Customers/1e6, color=clrs)\n",
    "axes[0,2].set_title(\"Total Number of Customers per Store Type (in Millions) \\n Fig 1.3\")\n",
    "\n",
    "axes[1,0].bar(df_train_store_mlc.groupby(by=\"StoreType\").sum().Customers.index,\n",
    "              df_train_store_mlc.groupby(by=\"StoreType\").Sales.mean(), color=clrs)\n",
    "axes[1,0].set_title(\"Average Sales per Store Type \\n Fig 1.4\")\n",
    "\n",
    "axes[1,1].bar(df_train_store_mlc.groupby(by=\"StoreType\").sum().Customers.index,\n",
    "              df_train_store_mlc.groupby(by=\"StoreType\").Customers.mean(), color=clrs)\n",
    "axes[1,1].set_title(\"Average Number of Customers per Store Type \\n Fig 1.5\")\n",
    "\n",
    "axes[1,2].bar(df_train_store_mlc.groupby(by=\"StoreType\").sum().Sales.index,\n",
    "              df_train_store_mlc.groupby(by=\"StoreType\").SalePerCustomer.mean(), color=clrs)\n",
    "axes[1,2].set_title(\"Average Spending per Customer in each Store Type \\n Fig 1.6\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From Fig 1.1 - Fig 1.3 We can see that Storetype A has the highest number of stores, sales and customers \n",
    "\n",
    "From Fig 1.4 - Fig 1.5 we got that Storetype B has the highest average Sales and highest average Number of Customers. \n",
    "\n",
    "Figure 1.6 tells us who has the highest average spending per Customer - it's StoreType D. Perhaps, this happened because no more similair shops around.\n",
    "\n",
    "Let's describe an **Assortment** level:\n",
    "\n",
    "> ### Analysis of Assortments with respect to Store Type :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.countplot(x=\"StoreType\",hue=\"Assortment\",order=[\"a\",\"b\",\"c\",\"d\"], \n",
    "                                     data=store_mlc,\n",
    "                                     palette=sns.color_palette(\"Set1\", n_colors=3))\\\n",
    "                          .set_title(\"Number of Different Assortments per Store Type\")\n",
    "\n",
    "store_mlc.groupby(by=[\"StoreType\",\"Assortment\"]).Assortment.count()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From both barplot and counting it is clearly to see. that almosr all of the stores have either 'a' assortment type or 'c', and 'b' type is only 9 stores!\n",
    "\n",
    "\n",
    "Talking about extra variety of products, we can say, that store type 'D' has the highest Sales per customer average actually and has mostly 'c' assortment type. This leads us to the reason of for having this high average in Sales per customer. Having variery in stores always increases the customers spending pattern.\n",
    "\n",
    "\n",
    "Also we can state another important fact: store type 'B' is the only one who has the 'b' assortment type (\"extra variety\"), looking at fig 1.4 and 1.5 exactly this store type has the highest number of customers and sales. It means, that having in the store not so much variety of products like 'c' assortment and not too basic like 'a' assortment will reach the high traffic in this store.\n",
    "\n",
    "\n",
    "We already saw the connection between **Sales and date in chapter 2**, there was common usuall pattern. Next step in exploring this theme is to look on connection between **Sales and Promotion withrespect to date**.\n",
    "\n",
    "> ### Analysis of Sales and Promotion withrespect to date:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"\\n\\t\\t\\t Sales with Promotion1,Promotion2 with respect to date\\n\")\n",
    "sns.factorplot(data = df_train_store_mlc, x =\"Month\", y = \"Sales\", col = 'Promo', hue = 'Promo2', row = \"Year\", \n",
    "               sharex=False, palette=sns.color_palette(\"Set1\", n_colors=2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"\\n\\t\\t\\t SalePerCustomer with Promotion1,Promotion2 with respect to date\\n\")\n",
    "sns.factorplot(data = df_train_store_mlc, x =\"Month\", y = \"SalePerCustomer\", col = 'Promo', hue = 'Promo2', \n",
    "               row = \"Year\", sharex=False, palette=sns.color_palette(\"Set1\", n_colors=2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is easy to see, here is a considerable difference between having promotion and not having. It means, that promotion really helps stores to change their amounts of sales to greater one.\n",
    "\n",
    "Another fact we can claim is the Promo2 variable, indicating a contunious promotion, works as \"inverse\" factor of promotion: in general, when there is no consecutive promotion, stores tend to sell more, then with consecutive promotion. Not so obvious, but we can suppose, that continous promotion is added for stores with very low sales to attract to attend customers. So those continuous promotion shows a tremending increase in the buying energy of customers.\n",
    "\n",
    "\n",
    "Looking through years there is a slight increase of total sales year over year. And as we already told on chapter2, there is a similair pattern in the months over the years with major pikes around the Christmas and Easter periods, in the 2/3 of Summer, and in March and April.\n",
    "\n",
    "And finally look at **Sales per day of week with respect to Promo**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.factorplot(data = df_train_store_mlc, x =\"DayOfWeek\", y = \"Sales\", hue='Promo' , sharex=False, \n",
    "               palette=sns.color_palette(\"Set1\", n_colors=2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (\"Number of Stores are opened on Sundays: \", \n",
    "       df_train_store_mlc[(df_train_store_mlc.Open == 1) & (df_train_store_mlc.DayOfWeek == 7)]\n",
    "       ['Store'].unique().shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No Promo on holidays. And on sundays work very few of stores. So this peak on 7th day of the week is normal, because of, perhaps, free time people having. So that's really important for stores to be opened on Sundays.\n",
    "\n",
    "Also there is a peak on Mondays with promotions, a tiny peak on Friday before the weekend."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ### Correlation Analysis \n",
    "\n",
    "First of all, we need numerical variables for our correlation analysis, so we have to transform not numerical variables to a numerical representation without leaking their logic.\n",
    "\n",
    "We remember that we have such variable as Open column. Since we have deleted all rows with zeros in this column (because of no profit for future model), now we can remove this column at all. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train_store_mlc.drop(columns='Open', inplace=True)\n",
    "\n",
    "# also in test : \n",
    "df_test_store_mlc.drop(columns='Open', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Let's check what types of variables we have : \\n\\n\", df_train_store_mlc.dtypes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have PromoInterval, StoreType, Assortment and StateHoliday as Obejcts.  \n",
    "\n",
    "For correlation functions, we need to make float variables. So we will create additional DF for this analysis.\n",
    "\n",
    "Thus we will transform it to category type, then we will get codes and set up type as float : \n",
    "\n",
    "But for modeling we will get dummies from categories - so just transform objects to categoires.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train_store_mlc['StateHoliday'] = df_train_store_mlc['StateHoliday'].astype('category')\n",
    "df_train_store_mlc['Assortment'] = df_train_store_mlc['Assortment'].astype('category')\n",
    "df_train_store_mlc['StoreType'] = df_train_store_mlc['StoreType'].astype('category')\n",
    "df_train_store_mlc['PromoInterval'] = df_train_store_mlc['PromoInterval'].astype('category')\n",
    "# there are some 0 are INT, so we will convert it to str : \n",
    "df_train_store_mlc.StateHoliday = df_train_store_mlc.StateHoliday.apply(lambda x: '0' if x==0 else x)\n",
    "\n",
    "\n",
    "df_correlation = df_train_store_mlc[['Store', 'DayOfWeek', 'Sales', 'Customers', 'Promo', 'SchoolHoliday', \n",
    "                               'CompetitionDistance', 'CompetitionOpenSinceMonth', 'CompetitionOpenSinceYear', \n",
    "                               'Promo2', 'Promo2SinceWeek', 'Promo2SinceYear', 'SalePerCustomer', 'Month', \n",
    "                               'Year', 'Day', 'StateHoliday', 'Assortment', 'StoreType', 'PromoInterval']]\n",
    "\n",
    "df_correlation['StateHoliday_cat'] = df_correlation['StateHoliday'].cat.codes.astype('float')\n",
    "df_correlation['Assortment_cat'] = df_correlation['Assortment'].cat.codes.astype('float')\n",
    "df_correlation['StoreType_cat'] = df_correlation['StoreType'].cat.codes.astype('float')\n",
    "df_correlation['PromoInterval_cat'] = df_correlation['PromoInterval'].cat.codes.astype('float')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the same for test : \n",
    "df_test_store_mlc['StateHoliday'] = df_test_store_mlc['StateHoliday'].astype('category')\n",
    "df_test_store_mlc['Assortment'] = df_test_store_mlc['Assortment'].astype('category')\n",
    "df_test_store_mlc['StoreType'] = df_test_store_mlc['StoreType'].astype('category')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "upper_triangle = np.zeros_like(df_correlation.corr(), dtype = np.bool)\n",
    "upper_triangle[np.triu_indices_from(upper_triangle)] = True \n",
    "\n",
    "f, ax = plt.subplots(figsize = (15, 10))\n",
    "sns.heatmap(df_correlation.corr(), ax=ax, mask=upper_triangle, annot=True, fmt='.2f', \n",
    "            linewidths=0.5, cmap=sns.diverging_palette(-100, 100, as_cmap=True))\n",
    "\n",
    "# cleaning\n",
    "del df_correlation\n",
    "gc.collect()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> So what have we got? \n",
    "\n",
    "There is a huge correlation with value 0.86 between Customers and Sales as we supposed before.\n",
    "\n",
    "Also significant value has Promo with Sales =0.4, since running a promotion increases that number .\n",
    "\n",
    "We can talk also about correlation between SalePerCustomers and Sales, but it's obvious and this variable is not available in test set.\n",
    "\n",
    "Furthermore, Assortment category with Sales has the 3rd place : 0.12\n",
    "\n",
    "Finally, we can see that StoreType does play a major role with SalePerCustomer = 0.43, this is probably due to the encoding of the store type variable.\n",
    "\n",
    "At this stage, we got a solid understanding of the distributions, the statistical properties and the relationships of our variables. The next step is to get clear with the metric for our model. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">>> ## Part 5. Metrics selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the most difficult issues is the question of choosing the most appropriate measures out of the groups. Due to the fact that each error measure has the disadvantages that can lead to inaccurate evaluation of the forecasting results, it is impossible to choose only one measure. \n",
    "\n",
    "In this competition we already have metric for model evaluation. Its' RMSPE. Let's consider it closer.\n",
    "\n",
    "In the begining there was definition: Root Mean Square Percentage Error (RMSPE). \n",
    "<center>\n",
    "<img src=\"https://habrastorage.org/webt/kq/pz/hx/kqpzhxo3app_0yojqxkxa42uzx0.png\" />\n",
    "where ${y_i}$ denotes the sales of a single store on a single day and ${y(hat)_i}$ denotes the corresponding prediction. Any day and store with 0 sales is ignored in scoring. \n",
    "\n",
    "\n",
    "We will take RMSPE, because of several reasons : \n",
    "On of the reasons is that this measure don't  depend on the magnitudes of the demand variables being predicted. \n",
    "\n",
    "Usually there is an appearance division by zero when the actual value is equal to zero. But as we remember, we already have got rid of all zeros in dataset. So it's not a problem for us now.\n",
    "\n",
    "So what about interpretation of this measure? \n",
    "\n",
    "> RMSPE: Forecasting Power\n",
    "\n",
    "> < 10%:\tHighly accurate forecasting\n",
    "\n",
    "> 10–20%:\tGood forecasting\n",
    "\n",
    "> 20–50%:\tResonable forecasting\n",
    "\n",
    "> more than 50%:\tInaccurate forecasting\n",
    "\n",
    "As you remember, I transformed the target variable Sales, SalePerCustomer, and Customers\n",
    ": they were coxboxed. So let's define function for inverse coxbox, using lambdas_coxbox defined early:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lambdas_coxbox: [0] Sales, [1] Customers, [2] SalePerCustomer\n",
    "\n",
    "\"\"\"\n",
    "# troubles with xgboost\n",
    "\n",
    "def ToWeight(y):\n",
    "    w = np.zeros(y.num_row(), dtype=float)\n",
    "    ind = y != 0\n",
    "    w[ind] = 1./(y[ind]**2)\n",
    "    return w\n",
    "\n",
    "\n",
    "def rmspe(yhat, y):\n",
    "    w = ToWeight(y)\n",
    "    rmspe = np.sqrt(np.mean( w * (y - yhat)**2 ))\n",
    "    return rmspe\n",
    "\"\"\"\n",
    "\n",
    "def rmspe(y, yhat):\n",
    "    return np.sqrt(np.mean((1 - yhat/y) ** 2))\n",
    "\n",
    "\n",
    "def rmspe_inv_cb(y, yhat, lambda_coxbox):\n",
    "    \n",
    "    y = inv_boxcox(y, lambda_coxbox)\n",
    "    yhat = inv_boxcox(yhat, lambda_coxbox)\n",
    "    \n",
    "    return rmspe(y, yhat)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">>> ## Part 6. Model selection\n",
    "\n",
    "For this competition we will try : RandomForestRegressor. \n",
    "\n",
    "This models should have good results because of the data: it's based on time series. \n",
    "\n",
    "So firstly, we will train models and fit using GridSearch for hyperparameter optimization to make a prediction. Also we will check results on hold-out sets. \n",
    "\n",
    "For sure, we can explore such models as LightGBM or NeuralNetworks, but for starting point it's enough.\n",
    "\n",
    "> RandomForest is a machine learning alogrithm used for classification and regression that is usually used with structured and tabular data. Its advantages:\n",
    "\n",
    "Random forest runtimes are quite fast, and they are able to deal with unbalanced and missing data.\n",
    "The process of averaging or combining the results of different decision trees helps to overcome the problem of overfitting.\n",
    "They also do not require preparation of the input data. You do not have to scale the data.\n",
    "\n",
    "> Its disadvantages:\n",
    "\n",
    "The main drawback of Random Forests is the model size. You could easily end up with a forest that takes hundreds of megabytes of memory and is slow to evaluate.\n",
    "They get a bit harder to interpret than regular deicison trees, since we are constructing of forest of more than 50 decision trees and more using grid search.\n",
    "\n",
    "\n",
    "P.S. I wanted to compare results with Xgboost, but some troubles appeared: it's hard to get values from xgboost Dmatrix - so no simple way to define user scoring function. So I left it for next time, and here I will explore only Random Forest. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# RandomForestRegressor\n",
    "rfr = RandomForestRegressor(n_estimators=10, \n",
    "                             criterion= 'mse', \n",
    "                             max_depth=5, \n",
    "                             n_jobs=7,\n",
    "                             random_state=17, \n",
    "                             warm_start=False)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">>> ## Part 7. Data preprocessing\n",
    "\n",
    "While exploring data and finding its percularities - we have already made some preprocessing : \n",
    "\n",
    "- Done work with all NANs: \n",
    " - - In additional store : 'CompetitionOpenSinceMonth', 'CompetitionOpenSinceYear' - Multilabel classification inplace of NANs.\n",
    " - - In test set Open NANs were replaced with 1 because of analysis particular store.\n",
    " - - Competition distance NANs were replaced with median of this variable. \n",
    " \n",
    "- Removed redudant column (Open), firstly got rid of zero rows in train set.\n",
    "\n",
    "- Coxboxed target Sales and really important variable Customers.\n",
    "\n",
    "- Some variables were transformed to categories.\n",
    "\n",
    "- Removed outliers. \n",
    "\n",
    "- Created new features such a Day, Month, etc.\n",
    "\n",
    "- Merged train and test sets with additional info in store DF. \n",
    "\n",
    "- All cleaning actions that we made for train set - we made for test set similarly.\n",
    "\n",
    "Now we will continue creating some new features and cleaning of the data.\n",
    "\n",
    "Firstly, let's condsider this columns : \n",
    "\n",
    "> Since the competition variables **CompetitionOpenSinceYear and CompeitionOpenSinceMonth** have the same underlying meaning, merging them into one variable that we call CompetitionOpenSince makes easier for the algorithm to understand the pattern and creates less branches and thus complex trees.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train_store_mlc['CompetitionOpenSince'] = np.where((df_train_store_mlc['CompetitionOpenSinceMonth']==0) & \n",
    "                                  (df_train_store_mlc['CompetitionOpenSinceYear']==0) , 0,\n",
    "                                  (df_train_store_mlc.Month - df_train_store_mlc.CompetitionOpenSinceMonth) + \n",
    "                                  (12 * (df_train_store_mlc.Year - df_train_store_mlc.CompetitionOpenSinceYear)))\n",
    "df_train_store_mlc.drop(columns=['CompetitionOpenSinceYear', 'CompetitionOpenSinceMonth'], inplace=True)\n",
    "\n",
    "# the same for test \n",
    "df_test_store_mlc['CompetitionOpenSince'] = np.where((df_test_store_mlc['CompetitionOpenSinceMonth']==0) & \n",
    "                                  (df_test_store_mlc['CompetitionOpenSinceYear']==0) , 0,\n",
    "                                  (df_test_store_mlc.Month - df_test_store_mlc.CompetitionOpenSinceMonth) + \n",
    "                                  (12 * (df_test_store_mlc.Year - df_test_store_mlc.CompetitionOpenSinceYear)))\n",
    "df_test_store_mlc.drop(columns=['CompetitionOpenSinceYear', 'CompetitionOpenSinceMonth'], inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **The StateHoliday** is not very important to distinguish (what type of holiday) and can be merged in a binary variable called is_holiday_state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train_store_mlc[\"is_holiday_state\"] = df_train_store_mlc['StateHoliday'].map({\"0\": 0, \"a\": 1, \"b\": 1, \"c\": 1})\n",
    "df_train_store_mlc.drop(columns=['StateHoliday'], inplace=True)\n",
    "\n",
    "# the same for test \n",
    "df_test_store_mlc[\"is_holiday_state\"] = df_test_store_mlc['StateHoliday'].map({\"0\": 0, \"a\": 1, \"b\": 1, \"c\": 1})\n",
    "df_test_store_mlc.drop(columns=['StateHoliday'], inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> It's better when working with decision tree based models to have dummy variables instead of categorical with different levels, because this alters the bias of the algorithm who will favor a higher weight to the categories like 4 and deprioritize levels like 1. And this problem could rise in the variables **Assortment and StoreType and PromoInterval**. So we will use the get_dummies function to instead do a binary encoding , as I said before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train_store_mlc = pd.get_dummies(df_train_store_mlc, columns=[\"Assortment\", \"StoreType\",\"PromoInterval\"], \n",
    "                              prefix=[\"is_Assortment\", \"is_StoreType\",\"is_PromoInteval\"])\n",
    "\n",
    "\n",
    "# the same for test \n",
    "df_test_store_mlc = pd.get_dummies(df_test_store_mlc, columns=[\"Assortment\", \"StoreType\",\"PromoInterval\"], \n",
    "                              prefix=[\"is_Assortment\", \"is_StoreType\",\"is_PromoInteval\"])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Also let's get dummies from **Day of Week**: \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_train_store_mlc  = pd.get_dummies(df_train_store_mlc, columns=['DayOfWeek'], prefix='Day')\n",
    "\n",
    "# the same for test \n",
    "df_test_store_mlc = pd.get_dummies(df_test_store_mlc, columns=['DayOfWeek'],prefix='Day')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# separating variables and target . dropping target from train .\n",
    "y = df_train_store_mlc.Sales.values\n",
    "all_vars_close_to_target = df_train_store_mlc[['Sales', 'Customers', 'SalePerCustomer']]\n",
    "\n",
    "df_train_store_mlc.drop(columns=['Sales', 'Customers', 'SalePerCustomer'], inplace=True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Now we will divide out train set on **80% training and 20% validation**.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_train_test, y_train, y_train_test = train_test_split(df_train_store_mlc, y, \n",
    "                                                                test_size=0.2, random_state=17)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">>> ## Part 8. Cross-validation and adjustment of model hyperparameters\n",
    "\n",
    "\n",
    "Let's start with RandomForestRegressor that we defined before.\n",
    "\n",
    "The\tmain parameters we will tune for RF is the number of trees and the size of the random subsets of features to consider when splitting a node. Also we will consider max_depth.\n",
    "\n",
    "We will use\t10 fold cross validation to get RMSPE while varying these parameters. Choosing 10 K-Folds makes sure we went through all of the data and didn't miss any pattern. \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we will check here what will be for initial RFR. Just for starting point. \n",
    "\n",
    "rfr = RandomForestRegressor(n_estimators=10, \n",
    "                             criterion='mse', \n",
    "                             max_depth=5, \n",
    "                             min_samples_split=2, \n",
    "                             min_samples_leaf=1, \n",
    "                             min_weight_fraction_leaf=0.0, \n",
    "                             max_features='auto', \n",
    "                             max_leaf_nodes=None, \n",
    "                             min_impurity_decrease=0.0, \n",
    "                             min_impurity_split=None, \n",
    "                             bootstrap=True, \n",
    "                             oob_score=False,\n",
    "                             n_jobs=4,\n",
    "                             random_state=31, \n",
    "                             verbose=0, \n",
    "                             warm_start=False)\n",
    "\n",
    "rfr_fin.fit(df_train_store_mlc, y)\n",
    "pred_1 = rfr_fin.predict(df_test_store_mlc)\n",
    "\n",
    "pred_2_inv_cb = inv_boxcox(pred_1, lambdas_coxbox[0])\n",
    "sample_submission_1 = sample_submission.copy()\n",
    "sample_submission_1.Sales = pred_1_inv_cb\n",
    "\n",
    "sample_submission_1.to_csv(\"subm1.csv\", index=False)\n",
    "\n",
    "# ok Here we got on LB : 0.36. Bad results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's tune! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = {'max_depth':(4,6,8,10,12,14,16,20),\n",
    "         'n_estimators':(4,8,16,24,48,72,96,128),\n",
    "         'min_samples_split':(2,4,6,8,10)}\n",
    "\n",
    "grid = RandomizedSearchCV(estimator=rfr,param_distributions=params,cv=10) \n",
    "grid.fit(X_train, y_train)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This is the best combination : \n",
    "\n",
    "grid.best_params_, grid.best_score_\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> We will check this model for valid set and test on LB in next chapters "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">>> ## Part 9. Creation of new features and description of this process\n",
    "\n",
    "We already created many new features before, here we will look on their importances.\n",
    "\n",
    "Let's consider feature importances:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(rfr_fin.feature_importances_)\n",
    "plt.figure(figsize=(18, 14))\n",
    "(pd.Series(rfr_fin.feature_importances_, index=X_train.columns)\n",
    "   .nlargest(20)\n",
    "   .plot(kind='barh'))\n",
    "plt.title('Feature Importance') # top 20 features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our top 5-6 most important variables are:\n",
    "\n",
    "1-Competitor Distance: This indeed impacts a lot the sales of a store like we saw previously in our EDA,when competition is very far stores tend to sell a lot more.\n",
    "\n",
    "2-Promo: Promotion is primordial for a store to increase its sales, it allows price breaking and thus more customers intersted in buying.\n",
    "\n",
    "3-Store: The Store itself represents a unique identificator for the algorithm to recognise which store has what attributes and indeed better accounts for the forecasting of those same stores in a future timeline.\n",
    "\n",
    "4-CompetitionOpenSince: The merging of this variable paid out and allowed us to give more accurate predicitions of the sales based on the time of opening of those competitors.\n",
    "\n",
    "5-6 -Promo2 Since week and Year: maybe it would better to convert them to one variable. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">>> ## Part 10. Plotting training and validation curves "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Got this functions from git and a bit improved\n",
    "\n",
    "\n",
    "def plot_with_err(x, data, **kwargs):\n",
    "    mu, std = data.mean(1), data.std(1)\n",
    "    lines = plt.plot(x, mu, '-', **kwargs)\n",
    "    plt.fill_between(x, mu - std, mu + std, edgecolor='none',\n",
    "                     facecolor=lines[0].get_color(), alpha=0.2)\n",
    "    \n",
    "def plot_validation_curve_gscv(gscv, param, logx=False):\n",
    "    plot_func = plt.semilogx if logx else plt.plot\n",
    "    res = gscv.cv_results_\n",
    "    x = [p[param] for p in res['params']]\n",
    "\n",
    "    # minus - neg-mse to mse\n",
    "    mu, std = -res['mean_train_score'], res['std_train_score']\n",
    "    lines = plot_func(x, mu, label='train')\n",
    "    plt.fill_between(x, mu - std, mu + std, edgecolor='none', facecolor=lines[0].get_color(), alpha=0.2)\n",
    "\n",
    "    # minus - neg-mse to mse\n",
    "    mu, std = -res['mean_test_score'], res['std_test_score']\n",
    "    lines = plot_func(x, mu, label='test')\n",
    "    plt.fill_between(x, mu - std, mu + std, edgecolor='none', facecolor=lines[0].get_color(), alpha=0.2)\n",
    "    plt.legend()\n",
    "    plt.grid(True)\n",
    "    \n",
    "def plot_learning_curve(reg, train_sizes, X, y):\n",
    "    scoring_fun = make_scorer(rmspe, greater_is_better=False)\n",
    "    N_train, val_train, val_test = learning_curve(reg,  X, y, \n",
    "                                                  train_sizes=train_sizes, \n",
    "                                                  cv=KFold(n_splits=3, shuffle=True, random_state=17),\n",
    "                                                  n_jobs=-1, scoring=scoring_fun)\n",
    "    \n",
    "    plot_with_err(N_train, -val_train, label='training scores')\n",
    "    plot_with_err(N_train, -val_test, label='validation scores')\n",
    "    plt.xlabel('Training Set Size'); \n",
    "    plt.ylabel('rmspe')\n",
    "    plt.legend()\n",
    "    plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_validation_curve_gscv(grid, 'max_depth')\n",
    "plot_validation_curve_gscv(grid, 'n_estimators')\n",
    "plot_validation_curve_gscv(grid, 'min_samples_split')\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_sizes = np.linspace(0.05, 1, 10)\n",
    "plot_learning_curve(rfr, train_sizes, X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">>> ## Part 11. Prediction for test or hold-out samples \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rfr_val=RandomForestRegressor(n_estimators=128, \n",
    "                             criterion='mse', \n",
    "                             max_depth=20, \n",
    "                             min_samples_split=10, \n",
    "                             min_samples_leaf=1, \n",
    "                             min_weight_fraction_leaf=0.0, \n",
    "                             max_features='auto', \n",
    "                             max_leaf_nodes=None, \n",
    "                             min_impurity_decrease=0.0, \n",
    "                             min_impurity_split=None, \n",
    "                             bootstrap=True, \n",
    "                             oob_score=False,\n",
    "                             n_jobs=7, \n",
    "                             random_state=17, \n",
    "                             verbose=0, \n",
    "                             warm_start=False)\n",
    "model_RF_test=rfr_val.fit(X_train,y_train)\n",
    "\n",
    "yhat=model_RF_test.predict(X_train_test)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "error=rmspe_inv_cb(y_train_test, yhat, lambdas_coxbox[0])\n",
    "print(round(error, 5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> #### Valid score : 0.17903"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The results on the LB score "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_test_store1 = df_test_store_mlc.copy()\n",
    "kaggle_yhat = model_RF_test.predict(df_test_store1)\n",
    "\n",
    "kaggle_preds = pd.DataFrame({'Id': sample_submission['Id'], 'Sales': inv_boxcox(kaggle_yhat, lambdas_coxbox[0])})\n",
    "kaggle_preds.to_csv(\"submission_3.csv\", index = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#kaggle_preds_ = pd.read_csv(\"submission_3.csv\")\n",
    "#kaggle_preds_\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "> #### Private Score : 0.16961\n",
    "\n",
    "> #### Public Score: 0.15653"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Nice results I think. But much more we can implement here. The best solution of this competition was around 0.1. So it encourages me. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">>> ## Part 12. Conclusions\n",
    "\n",
    "\n",
    "We made pretty huge data analysis. \n",
    "\n",
    "Some interesting ideas were implemented such as multilabel classification for NANs, and here we have really big space for improving model.\n",
    "\n",
    "Unfortunately, I didn't have enough time to deal with some problems with Xgboost, so it's good way developing and improving the solution. So we explored only RandomForest, and next points will be : discover results of Xgboost, LightGBN and perhaps NN. So it's really wotrh trying more models next time. And then compare the results. \n",
    "\n",
    "So maybe I will create some kind of own Xgboost and then test the results :)\n",
    "\n",
    "We can understand from this project the flexibility and robustness of a decision tree based model like RandomForest which helped us predict the Store Sales of Rossman based on attributes that defines each store and its surroundings.\n",
    "\n",
    "As we can see, it always delivers a good predicition score while not having a lot of modifications and difficulties capturing the patterns hidden in the data. Fortunately we had a train set that was large enough for it to converge but in general RandomForest performs not so bad on small sets since its resampling method (bagging) and its random production of trees allow the bias to remain not so high. \n",
    "\n",
    "\n",
    "I can notice that there are many ways to build new features - and it's really interesting and challenging. This steps that I covered - are the basis for my future analyses. And also I think that using hyperparameter optimization techniques like Gridsearch and RandomizedSearch is crucial to any Machine Learning problem since it allows the algorithm to not just limit itself on its defaulted parameters but to discover new opportunities of combining those parameters to reach a local optima while training on the data.\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "project_template.ipynb",
   "provenance": [],
   "version": "0.3.2"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
